### LOAD DATA

# A) Import Results

# Package to read .mat files:

library(R.matlab)

# Initialize storage (parameters and obj. fn. value):

param_obj <- numeric(0)

# Load first file:

data <- readMat("mh_search01.mat")

# Add to storage:

param_obj <- rbind(param_obj, cbind(data$parameters, data$densities) )

# Load another file:

data <- readMat("mh_search02.mat")

# Add to storage:

param_obj <- rbind(param_obj, cbind(data$parameters, data$densities) )

# B) Recovered Set ("equilibrium points")

# Log pseudo-density is 0:

eq.points <- param_obj[param_obj[,5] == 0, ]  

# Simplify to unique points and format object:

eq.points <- unique(eq.points)
eq.points <- as.data.frame(eq.points)

# Variable names:

colnames(eq.points) <- c("b","c0","cN","cT","obj")

# Save:
# (uncomment one file name below, as appropriate)

rm(data)

# save.image("mh_search_c2_pop.RData")
# save.image("mh_search_c2_samp.RData")



### PLOT RESULTS

# (after all MCMC chains have been imported and both files above have been saved)

# Reload population results:

load("mh_search_c2_pop.RData")

# Reload sample results:

samp.50 <- new.env()

load("mh_search_c2_samp.RData", envir=samp.50)

###

# A) All Values of Beta (Figure A2)

# Plot b-cN:

png("b-cN_dif.png", width=360, height=360)

attach(samp.50$eq.points)
plot(b,cN, pch=20, xlim=c(-2,+2), ylim=c(-2,+2) , col="gray")
detach(samp.50$eq.points)

attach(eq.points)
points(b,cN, pch=20, xlim=c(-2,+2), ylim=c(-2,+2) )
detach(eq.points)

dev.off()

# Plot b-cT:

png("b-cT_dif.png", width=360, height=360)

attach(samp.50$eq.points)
plot(b,cT, pch=20, xlim=c(-2,+2), ylim=c(-1,+3) , col="gray")
detach(samp.50$eq.points)

attach(eq.points)
points(b,cT, pch=20, xlim=c(-2,+2), ylim=c(-1,+3) )
detach(eq.points)

dev.off()

### 

# B) Positive Beta (Figure 3)

# Parameter restriction:

attach(eq.points)
eq.points.pos <- eq.points[b >= 0, ]
detach(eq.points)

attach(samp.50$eq.points)
samp.50$eq.points.pos <- samp.50$eq.points[b >= 0, ]
detach(samp.50$eq.points)

# Table 4:

round(cbind(
	apply(eq.points.pos, 2, min), apply(eq.points.pos, 2, max),
	apply(samp.50$eq.points.pos, 2, min), apply(samp.50$eq.points.pos, 2, max)
), 2)

# Plot b-cN:

png("b-cN_pos_dif.png", width=360, height=360)

attach(samp.50$eq.points.pos)
plot(b,cN, pch=20, xlim=c(-0.5,+1.5), ylim=c(-1.5,+1.0), col="gray")
detach(samp.50$eq.points.pos)

attach(eq.points.pos)
points(b,cN, pch=20, xlim=c(-0.5,+1.5), ylim=c(-1.5,+1.0) )
detach(eq.points.pos)

dev.off()

# Plot b-cT:

png("b-cT_pos_dif.png", width=360, height=360)

attach(samp.50$eq.points.pos)
plot(b,cT, pch=20, xlim=c(-0.5,+1.5), ylim=c(-1.0,+1.5), col="gray")
detach(samp.50$eq.points.pos)

attach(eq.points.pos)
points(b,cT, pch=20, xlim=c(-0.5,+1.5), ylim=c(-1.0,+1.5) )
detach(eq.points.pos)

dev.off()
